#### PREAMBLE ####

library(ggplot2)
library(geosphere)
library(dplyr)
library(data.table)
library(tidyr)
library(zoo)
library(sf)
library(ggspatial)
library(RColorBrewer)
library(estimatr)
library(broom)
library(texreg)
library(did)


#### TEMPLATE ####
goldenScatterCAtheme <- theme(
  panel.background = element_rect(fill = "white"),
  aspect.ratio = ((1 + sqrt(5))/2)^(-1),
  axis.ticks.length = unit(0.5, "char"),
  axis.line.x.top = element_line(size = 0.2),
  axis.line.x.bottom = element_line(size = 0.2),
  axis.ticks.x = element_line(size = 0.2),
  axis.text.x = element_text(color = "black", size = 10),
  axis.title.x = element_text(size = 10, 
                              margin = margin(t = 7.5, r = 0, b = 0, l = 0)),
  axis.ticks.y = element_blank(),
  axis.text.y = element_text(color = "black", size = 10,
                             margin = margin(t = 0, r = -4, b = 0, l = 0)),
  axis.title.y = element_text(size = 10,
                              margin = margin(t = 0, r = 7.5, b = 0, l = 0)),
  #legend.key = element_rect(fill = NA, color = NA),
  panel.grid.major.x = element_blank(),
  panel.grid.major.y = element_line(color = "gray45", size = 0.2),
  strip.background = element_blank(),
  strip.text.x = element_text(size = 8),
  strip.text.y = element_blank(),
  strip.placement = "outside",
  panel.spacing.x = unit(1.25, "lines"),
  panel.spacing.y = unit(1, "lines")
)
class(goldenScatterCAtheme)
#### LOAD DATA ####

set.seed(20230711)

#load main data
df <- read.csv("Outputs/Analysis all muns.csv", stringsAsFactors = F, na.strings = c("", "NA"))
df$INEGI <- sprintf("%05d", df$INEGI)
df$INEGI <- as.character(df$INEGI)
df$time <- as.Date(df$time)

#### TRANSFORMATION ####

#states treated by US trade restrictions
unique(df$State[df$ustrade==1])
df$ustrade[df$time>=as.Date("2016-06-01") & df$trade==1 & (df$State=="MICHOACAN")] <- 1
#updating ustrade.
df$ustrade <- df$ustrade * df$trade

#### SUBSET TO PROPENSITY MATCHED SAMPLE ####

matched <- read.csv("Outputs/propensity_matched_muns.csv",stringsAsFactors = F)
matched$INEGI <- sprintf("%05d", matched$INEGI)
matched$INEGI <- as.character(matched$INEGI)

df <- df[df$INEGI %in% matched$INEGI,]


#### MAP OF TREAT_GROUP ####

library(ggplot2)
library(ggspatial)
library(ggpattern)

#create data of treat_group by municipality
map_data <- df[,c("INEGI", "treat_group")]
map_data <- map_data %>% 
  distinct(INEGI, .keep_all=T)

map_data$treat_group[map_data$treat_group=="Turnson"] <- "Export certified between 2011-2019"
map_data$treat_group[map_data$treat_group=="Alwayson"] <- "Export certified before 2011"

#municipalities shape file
mun <- st_read("Source data/Map files/00mun.shp")
st_crs(mun)
#states shape file
state <- st_read("Source data/Map files/Entidades_2010_5.shp")

#MERGE QUERY AND MUN DATASETS
mun <- left_join(mun, map_data, by = c("CVEGEO"="INEGI"))

mun$treat_group[is.na(mun$treat_group)] <- "Missing data or excluded"

#relevel treat group factor
map_data$treat_group <- factor(map_data$treat_group, levels=c("Missing data", "Control", "Export certified before 2011", "Export certified between 2011-2020"))

#pick colors
display.brewer.all(n=NULL, type="all", select=NULL, exact.n=TRUE)
colors <- brewer.pal(n=9, name="Pastel1")
white <- "white"
grey <- "gray85"
red <- colors[1]
green <- colors[3]
darkred <- brewer.pal(n=9, name="Reds")[5]
darkgreen <- brewer.pal(n=9, name="Greens")[5]
black <- "black"

#add centroids for plotting state names
state <- cbind(state, st_coordinates(st_centroid(state)))
#rename some states
state$NOM_ENT <- as.character(state$NOM_ENT)
state$NOM_ENT
state$NOM_ENT[3] <- "Mexico"
state$NOM_ENT[14] <- "Michoacan"
state$NOM_ENT[21] <- "Queretaro"
state$NOM_ENT[22] <- "San Luis Potosi"
state$NOM_ENT[27] <- "Nuevo Leon"
state$NOM_ENT[30] <- "Yucatan"
state$NOM_ENT[26] <- "Veracruz"
state$NOM_ENT[28] <- "Coahuila"
state$NOM_ENT[1] <- ""

mun <- st_transform(mun, crs=6372)
st_crs(mun)
state <- st_transform(state, crs=6372)
st_crs(state)

us_trade_muns <- unique(df$INEGI[df$ustrade==1])
mun$ustrade <- NA_character_
mun$ustrade <- ""
mun$ustrade[mun$CVEGEO %in% us_trade_muns] <- ", Able to export avocados to the U.S."
mun$treat_group <- paste0(mun$treat_group, mun$ustrade)

h <- 6
w <- 1.618*h
#create pdf of map with states and municipalities and query presence in 2010
#Plot the data
ggplot() +
  ggspatial::annotation_spatial(mun) +
  ggspatial::layer_spatial(mun, aes(fill=treat_group), color="gray14", linewidth=.05)+
  scale_fill_manual(values=c(grey, red, darkred, green, darkgreen, white))+
  geom_sf(data = state, fill=NA, color=black, linewidth=.4) +
  geom_text(data = state, aes(X, Y, label = NOM_ENT), size = 4)+
  theme_void()+
  coord_sf(xlim = c(2000000, 3000000), ylim = c(500000, 1200000), expand = FALSE)+
  theme(legend.position = c(0.2, 0.2), legend.title=element_blank(), legend.key=element_blank())
ggsave("Outputs/plots/Fig13_Mexico_state_mun_treat_groups_propensity_match.png", width=w, height=h)

#
